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We consider the problem of estimating spatially varying coefficients of 
partial differential equations from observation of the solution and of the 
right hand side of the equation. We assume that the observations are 
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estimation of coefficients in an elliptic and a parabolic partial differential 
equation. 
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1. Introduction 


This paper discusses the problem of estimating spatially varying 
coefficients of differential equations using observations of the solution and 
the right hand side of the equation. We are restricting the discussion to 
distributed observations. 

A common estimation procedure is to try to find parameters such that by 
solving the equation, the solution is as close as possible (in some norm) to 
the observations [1] ("output least squares"). Another approach is to look at 

the equation as an equation for the parameters [2]. The approach we are 

discussing here is closer to the second one: we seek a solution (for the 

coefficients) in the least squares sense. This approach and the first 
approach are two special cases of another approach which estimates the 

parameters by minimizing a weighted sum of the error between observed and 

computed solutions with the residual of the equation. That is, given the 

equation 

L (g)V = f 

one identifies g by minimizing over (g, V)C A the functional 

J (V,g) = II V - V*n 2 + a IIL (g)V - f*ll 2 

where V ,f are the observed quantities, and A is a set of admissible pairs of 
functions. 

Let (v(a), g(a)) be the minimum of J. Then o °°, (v(a), g(a)) 

converge to the solution of the "output least squares" solution. In the 

case a 0, (^(a), g(a)) converge to the solution of the minimization 
problem discussed in this paper ("equation error approach"). 
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The disadvantage of using the "equation error" approach is that one can 
treat only problems in which measurements are distributed in the domain and in 
which enough measurements are available. Also, this case may be more 
sensitive to noisy data than the "output least squares" approach because one 
applies L(g) to V , and if L(g) is a differential operator, a large error may 
be introduced even for a small amount of noise in V . 

Section 2 describes the mathematical formulation of the problem and some 
basic rules for discretization. A multigrid procedure for solving the 
resulting equations is presented in Section 3. In Section 4, we present some 
numerical computations in which we estimate coefficients of an elliptic and 
parabolic equation. The elliptic case appears in problems involving large 
space structures [1], and the parabolic case arises in oil exploration and 
recovery [1], 


2 . Formulation of the Problem and Discretization Method 

Let L(g) be a differential operator depending on a set of coefficients 
g(x) = ( gl (x),..., gjl (x)). We wish to estimate jg from the observed solution 
V(x) and a right hand side of the equation 

(2.1) L(£) V(x) = f(x) x e SI. 

The problem is generally ill-posed. In order that g_(x) will be 

identifiable, V(x) must satisfy some necessary conditions. (See [2] for 

example.) We assume throughout the paper that the observations V(x), f(x) 


are such that: 
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(Pl) observations are distributed in ft 

(P2) V(x), f(jc) are such that g(x) is identifiable 

(P3) g(x) depend continuously onV(x), f(x). 

In this paper, we are not discussing conditions on L(^), V(x), f(x) that 
imply (P2) and (P3). Basically, these assumptions imply that the 
identification problem is well posed. In a companion paper [5] , we discuss 
such conditions for the example given in section 4.1. Our approach for 
estimating j[(x) will be to look for J[(x) such that equation (2.1) will be 
satisfied in a least squares sense, i.e., 

(2,2) mln 11 L ^ V * “ f * 11 2 

g(x)€A 

where V , f are given in at discrete points, and A is a class of admissible 

parameters. The first step is to construct an interpolant V that is 

defined everywhere in ft. This will define a problem on the continuous 

level. Next we consider the question of discretization. 

ti li Vi 

Let the discrete domain be ft , on which V* , f* are given (by 

interpolation from the measurements). Let A* 1 , A^ be spaces of discrete 

functions defined on the discrete domains ft* 1 , ft** respectively, where 
H h H 

ft eft. A will be the set of admissible parameters. The dimension of 
A** may not be the same as that of A* 1 . The discretization of (2.1) is then 


L h (g H ) v£ 



g eA , 


*» 


f^eA h 


(2.3) 
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where L^(g^) is a matrix approximation to L(jf) . The identification 
procedure then takes the form 


(2.4) 


min H L^(g^) vj 1 

H.H 

g eA 


f h II 2 

* "2 


The method one uses for discretization is a crucial one for the quality of the 
estimated quantities. We suggest the following rule: 


H 

Discretization Rule : The number of parameters in j[ should be less than the 

number of equations in (2.3). 


The above rule is in accordance with results by others (for example [6]). 

This rule seems to guarantee that no spurious oscillations are developed in 
H H 

g . In practice, we take g to be defined on a grid twice as coarse as the 
ti h 

grid on which V* , f* are defined. 

I- 

Let the matrix B(V n ) be defined by 


(2.5) 


B(V h ) g H = L h (g H ) V h 


Note that, although L^(g H ) is a square matrix, B(V^) is a rectangular one 
H h 

since involves less parameters than V . With the above definition of 
B(V n ) , our minimization problem takes the form 


( 2 . 6 ) 


min 

H cA H 

g eA 


11 B 


(V* h )g H - 


This leads, in the case of linear L, to the following system of equations 
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(2.7) [B (V* h )] T B (V* h ) g H = [B(V^ h ) ] T f* h . 

We next describe an efficient algorithm for solving the discrete system 

(2.7) . 


3. Multigrid Procedure 

The main element of a good multigrid procedure is a relaxation scheme. 
We start by describing this part. 

Relaxation . Since [B (V^)] 1 B (V^ n ) is symmetric positive definite, it 
seems that Gauss-Seidel relaxation will be appropriate. However, for this we 
have to compute [B (V^) ] T B(V* h ) . For general problems where g H , V* h lie 
on different grids, this may be too involved. As an alternative, we suggest 
Kaczmarz relaxation for 

(3.1) B (V* h ) g H = f* h . 


From a theorem of Tanabe [3], this relaxation converges to the least squares 
solution plus the projection of the initial approximation onto the kernel of 

L iji 

[B (V* )] . Considering the task of programming the algorithm, it is much 

h X ti 

simpler to take this route over the one of computing [B (V* 1 )] 1 B (V*). 

The i-th step in Kaczmarz relaxation is given by: 

For j = 1 , . . . , n do: 


H 

g j 


H . , , 
g. + <5 . b. . 
3 1 iJ 


where 


B(v*) = (b 


. . ) and 6 . = 
ij i 


(£ - B(V* h )g H ) 


§ 'hj 
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We now come to the description of a multigrid cycle. We begin by 
describing a two-grid cycle. 

Two-grid Cycle 

v oh Vi 

Suppose we are given two grids ft, ft z where ft z C ft. Let 

L OL r or 

l* ^i (i = l>2) be spaces of grid functions defined on ft, ft z 

respectively. Assume the ft^-grid equations are 

B h g h = f h g h €^, f h €A^. 

Let be an interpolation operator from A 2 * 1 to A^, and 1^ = (l!j ) ’ 

/v# 

Given an approximate solution g to the above equation, a two-grid 
cycle for improving it consists of performing (A) - (D): 

(A) Relax the equation B^gb = f* 1 times, starting with 

g h yielding g h . 


(B) Solve (approximately) the coarse grid equation 


_2h 2h T 2h ,,h n h-h. 
Bg - I h (f - B g ); 


,2h 


- r 2h R h T h 

' h B *2h’ 


denote by 



the (approximate) solution. 


— h 
g «■ 


-h 

g 



~2h 

g • 


(C) Perform the correction 
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(D) Relax the equation B^g* 1 = f^ 1 ^ times, starting with g* 1 

yielding g . 

Note that is an operator from A^ to A^ and the definition given 

for it is only one possibility. 

In order to obtain full efficiency of such an algorithm one solves the 
coarse grid equations in step (B) by a similar algroithm using a still coarser 
grid. Applying this idea recursively we get the basic multigrid cycle which 
is defined next. 

Multigrid Cycle 

Given a sequence of discretizations with mesh sizes h,> h„ > ...h , on 

12m 

Ir If- 1 ]f 

grids Q , where h^ = 2h^ + ^ and H c & • Let the h^-grid equation 
be 

(3.3) B k g k = f k ; f k eA k , g k eA k , 

k k+1 k 

where B approximates B (k < m) and A^ (i=l,2) are spaces of 

|q 

functions defined on the h^-grid. Let an Interpolation operator 

a k- 1 k k- 1 k k- 1 

from A^ to A^, and 1^ be a restriction operator from A^ to A^ 

j /N^k 

The algorithm for improving a given approximate solution g to (3.3) is 
denoted by 

(3.4) £ k <- MG (k, g k , f k ) 


and is defined recursively as follows: 
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(1) If k = 1, solve (3.3) by several relaxations otherwise do steps (A) - 
(D). 


(A) Perform Vj relaxation sweeps on (3.3), resulting in a new approximation 

-k 

g • 


(B) Starting with g k ^=0, perform y 


times the following cycle 


~k-l 

g 
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A/](- 1 

MG(k-l,g R 



(f k ' B k g k )), 


(C) Calculate 


=k <• 
g 


-k 

g 


+ I 


k-1 


~k-l 

g 


(D) 


= lc 

Perform additional relaxation steps on (3.3) starting with g and 

yielding the final g k of (3.4). 


The cycle with y = 1 is called V(v^ ,v 2 ) -c ycle and the one with y = 2 
is called M^^v^). 

Full Multigrid Algorithm (FMG) 

In order to obtain full efficiency, it is better to start from a good 
initial guess. This can be obtained by solving a similar problem on a coarser 
level. Applying the idea recursively, we get the following algorithm which is 
called full multigrid algorithm and is denoted by N-FMG. 

1. Solve (3.3) for k=l using several relaxations. 
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2. k = k+1 

~k k ~k-l k 

g = ir^-i S , ^ iss an interpolation operator. 

3. Perform the cycle 

~k ,, ~k ,k> 

g «- MG (k, g , f ) 

N times. 

4. If k < m, go to 2; otherwise, stop. 

4. Numerical Examples 

In this section, we demonstrate the effectiveness of the method 
described. Results are given for two identification problems. In all cases 
the grids are finite and the spaces defined them are of grid functions (finite 
dimensional), with the 

4 . 1 Elliptic case 


Let g = g (x,y) ; (x,y) e = [0,1] x [0,1], and 

T , N _ 3 3 . 3 3 

L = 3x g 3x + 3y g 3y * 

Determine g(x,y) from observations of V(x,y), f(x,y) which satisfy 


L(g) V = f in ft 
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Discretization V, f were discretized on a uniform mesh 


0 = {(ih, jh), 0 < i, j < N} N • h = 1, 


g was discretized on a mesh fi , with H = 2h. The discrete equation is 


g H (JL+Li LVi) _ H (Vr!i±i) 

g i + V 2 > j h 2 ; g ±- l /2,3 u2 ; 


+ g H . ~ U ij } _ (VVj-l) - f 

g i. J+V h 2 } g i,j-V 2 ( .2 ) f ij 


where j, . and g 1 ? j, are given from the nodal values of g H on the 

'2 » J '2 

grid fi by linear interpolation. This gives approximately 4 times as many 
equations as there are unknowns. At every level of discretization, we 
maintained the relation between g H , U* 1 , i.e., g* 1 is always defined on a grid 
twice as coarse as the one for U h . The £ 2 norm of a function U h 
defined on a grid ft* 1 is given by 


" Uh,l 2 = h h i U Ji| 2 * 

(i,j)eo n 1J 


All interpolation was bilinear and residuals were transferred by the 9- 
point averaging operator. 


k-1 = J_ 
k 16 


1 

2 

LI 


2 

4 

2 
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The following tables summarize the numerical experiments. In all the examples 
1-FMG using W(2,l) cycle is shown. The numbers are taken at the end of the 1- 
FMG algorithm. In all the examples, the identification problem is well-posed. 


Example 1 V(x,y) = x(l-x) y(l-y) 

g(x,y) = 1 + x + y 


Level # 

h * 

II g " g « 2 

II residuals II ^ 

1 

.106 (-4) 

.370 (-5) 

2 

.328 (-6) 

.190 (-6) 

3 

.324 (-7) 

.409 (-7) 

4 

.0 

.0 (-9) 


Example 2 V(x,y) = sin(wx) sin(iry) 

g(x,y) = 1 + x + y 


Level // 

. h 
II g 

* 

- g l 2 

II residuals II ^ 

1 

.114 


.452 

(-4) 

2 

.289 

(-1) 

.322 

(-1) 

3 

.694 

(-2) 

.818 

(-2) 

4 

.178 

(-2) 

.242 

(-2) 
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Example 3 


V(x,y) = sin(irx) sln(iry) 


g(x,y) = 


1 + x 
1.5 


x < 

X > 


5 

5 


Example 


Level 

# 

h * 

n g - g « 2 

K residuals II 

2 

1 


. 8 1 2 ( — 1 ) 

.289(-5) 


2 


• 1 96( — 1 ) 

.439(-l) 


3 


.530(-2) 

.181(-1) 


4 


.153(-2) 

.603(-2) 



V(x,y) 

= x(l-x)y(l-y) 



g(x,y) 

) 1 + x, x < .5 

(1.5, x > .5 


Level 

# 

h * 

n g - h n 

n residuals II 

2 

1 


.206(-l ) 

.237 ( — 5 ) 


2 


.536(-2) 

.256(-2) 


3 


.204(-2) 

•959(-3) 


4 


.750(-3) 

. 3 1 4( — 4 ) 

_ 


4.2 Parabolic Case 

2 

Let g = (g(x,y), s(x,y) ); (x,y) e = [0,1] ; and 


Ug.s) - -s (X,y) | T + |_g|_ + | ?g |_ 
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Determine s(x,y), g(x,y) from observation of V(x,y,t) f(x,y,t) such that the 
following equation holds 

L (g,s) V = f (x, y) e ft , 0 < t < T . 

o 

In this problem, we used observed data at a few different times. We 

assumed that both V and V are given together with f. 

h M 

Let ft be defined as in the previous section and r be defined as 



2 

The total number of equations is n times the numbers of time observa- 
tions. In our example, five time observations were used. 

In defining coarse grids, we maintained the relation between v,f,s,g as 
on the fine level, i.e., v,f,s use a grid twice as fine as g uses. We 
coarsened in space only, leaving the number of time measurements to be the 
same on all levels. Intergrid transfers were exactly as in the elliptic case 


(sec 4.1). The relaxation consists of two steps. The first one was a 

H h 

Kaczmarz relaxation in which both g and s were relaxed. The second one was 
a pointwise relaxation of s n only. In this relaxation, the points 


(i,j) e ft were scanned in lexicographic ordering, and at each point 



was 
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changed to minimize 


„ •» h / H h. 
BL (g ,s ) 


v h * - f h n . 


The following tables summarize the numerical results. Results are given 
for 3-FMG using w(2,2) - cycle. 


Example 1 


v(x,y,t) = sin(Trx)sin(Try) (t+1) 
g(x,y) = 1 + x + y 
s(x,y) = 1. 


Level # 

Cycle # 

h * 
Bg "g 


h 

II s -s 

* 

II 

II residuals II 2 

1 

10 

.114 


8.11 

(-5) 

3.52 

(-5) 

2 

1 

2.75 

(-2) 

1.17 

(-2) 

4.64 

(-3) 


2 

2.75 

(-2) 

5.06 

(-3) 

4.14 

(-3) 


3 

2.75 

(-2) 

2.72 

(-3) 

3.81 

(-3) 


4 

2.75 

(-2) 

1.98 

(-3) 

3.54 

(-3) 

3 

1 

6.81 

(-3) 

1.96 

(-2) 

1.50 

(-3) 


2 

6.79 

(-3) 

8.75 

(-3) 

1.38 

(-3) 


3 

6.78 

(-3) 

3.95 

(-3) 

1.29 

(-3) 


4 

6.78 

(-3) 

1.95 

(-3) 

1.20 

(-3) 

L_ 

4 

1 

1.77 

(-3) 

3.30 

(-2) 

2.41 

(-3) 


2 

1.72 

(-3) 

1.98 

(-2) 

2.58 

(-3) 


3 

1.70 

(-3) 

1.20 

(-2) 

2.67 

(-3) 


4 

1.69 

(-3) 

7.96 

(-3) 

2.56 

(-3) 
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Example 2 


V(x,y,t) = sin(irx) sin(7ry)(t+l ) , s(x,y) = 1 



g(x,y) = 

1+x j x < .5 
1.5 |x > .5 




Level // 

Cycle // 

h * 
ng -g « 

„ h 
Its -s 

* 

# 

n residuals II ^ 

1 

10 

8.19 (-2) 

4.49 

(-5) 

1.93 (-5) 


| 

1 

1.90 (-2) 

6.98 

(-3) 

.137 


2 

1.90 (-2) 

1.01 

(-2) 

.128 


3 

1.91 (-2) 

1.46 

(-2) 

.119 


4 

1.92 (-2) 

1.88 

(-2) 

.111 


3 

1 

5.08 (-3) 

1.02 

(-2) 

5.78 (-2) 


2 

5.17 (-3) 

7.02 

(-3) 

5.37 (-2) 


3 

5.22 (-3) 

6.86 

(-3) 

4.99 (-2) 


4 

5.25 (-3) 

7.79 

(-3) 

4.66 (-2) 


4 

1 

1.43 (-3) 

1.68 

(-2) 

2.29 (-2) 


2 

1.44 (-3) 

1.01 

(-2) 

2.17 (-2) 


3 

1.45 (-3) 

6.83 

(-3) 

2.06 (-2) 


4 

1.46 (-3) 

5.59 

(-3) 

1.97 (-2) 


I 

4.3 Discussion 

These numerical examples clearly demonstrate the effectiveness of the 
method of discretization as well as the solution process. 
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2 

In the elliptic problem, the coefficients converge in a rate of 0(h ) 
in the smooth case. When derivatives of the coefficients are not smooth, a 
slower rate is obtained as one would expect. The multigrid algorithm solves 
the problem up to the level of discretization errors in just one multigrid 
cycle of the type FMG-W(2,1). 

2 

In the parabolic case, the coefficient g(x,y) shows 0(h ) convergence 
(in the smooth case), and this is obtained after the first cycle. s(x,y) 
behaves worse taking more cycles to reach the level of discretization error. 
This can be explained as follows: a change of order 0(h ) in g(x,y) may lead 
to a change of 0(h) in s(x,y). Hence s(x,y) may reach convergence only after 
g(x,y) has converged. It is possible that, although s(x,y) has not reached 
convergence, g(x,y) will be accurate up to discretization errors if the error 
in s(x,y) times V t (x,y,t) is of the level of truncation errors. The 
behavior of the residuals in examples of Section 4.2 is not clear. It goes up 
when going from level three to level four, and it seems that it should have 


gone down 
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